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^This  phase  of  the  study  concentrated  on  refinement  of  vegetation  extrac¬ 
tion  techniques  developed  in  phases  1  through  3.  Software  was  also  generated 
for  scaling  and  mosaicking  vegetation  data  derived  from  separate  adjacent 
scenes,  as  appropriate  for  covering  a  given  map  area.  The  raster  format  classi¬ 
fication  files  were  converted  to  polygon  vector  files  in  the  Standard  Interface 
Format. 


SUMMARY 


The  overall  Terrain  Data  Extraction  Study  is  concerned  with  the  development  of 
techniques  for  the  extraction  of  terrain  features  from  digital  aerial  imagery. 
Typical  features  in  the  vegetation  category  are  vegetation  boundaries  and  tree 
stem  spacing  in  forested  areas.  Extracted  data  pertaining  to  features  like 
these  are  an  input  to  decision-making  processes  associated  with  military  opera¬ 
tions  such  as  cross-country  movement.  Phase  4  effort,  reported  here,  has  built 
on  the  base  of  vegetation  feature  extraction  techniques  developed  in  the  previous 
phases. 

Automated  delineation  (boundary  outlines),  in  a  selected  digitized  aerial  panchro¬ 
matic  photograph,  has  been  performed  for  vegetation/landcover  classes,  intervals 
of  percent  canopy  closure,  and  intervals  of  tree  stem  spacing.  Images  have  been 
produced  showing  the  land-cover  class  and  interval  themes,  and  the  theme  boundaries 
have  been  extracted  and  recorded  on  magnetic  tape  in  the  Standard  Interface  Format. 

To  have  utility,  terrain  data  boundaries  must  be  produced  for  specified  map  areas 
normally  encompassing  multiple  frames  of  adjacent  aerial  imagery.  Algorithms 
have  been  developed,  therefore,  for  digital  scaling  and  mosaicking  of  terrain  data 
derived  from  separate,  adjacent  scenes.  Software  to  implement  the  algorithms  has 
also  been  produced. 

Within  a  mosaic  of  multiple  aerial  images,  the  terrain  data  boundaries  for  any 
given  category  are  likely  to  be  quite  complex.  Multiple  polygons,  polygons  with 
holes,  polygons  within  polygons,  etc.  are  involved,  with  each  polygon  containing 
a  large  number  of  vertices.  Algorithms  have  been  developed  to  automatically  and 
systematically  examine  the  terrain  data  binary  themes  and  to  produce  the  appro¬ 
priate  single-pixel-width  outlines  (polygons).  Procedures  also  have  been  developed 
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to  convert  the  outline  data  to  the  Standard  Interface  Format  suitable  for  input 
to  the  ETL  Interactive  Graphic  Design  System.  These  procedures  and  algorithms 
have  been  implemented  with  appropriate  software. 


A  detailed  design  for  an  interactive  system  capable  of  performing  the  above  tasks 
was  prepared  early  in  Phase  4  and  documented  in  a  separate  report  entitled  Inter¬ 
active  Image  Analysis  System  Design,  dated  and  submitted  to  USAETL  January  1983. 
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INTERACTIVE  DIGITAL  INAGE  PROCESSING  FOR  TERRAIN 


DATA  EXTRACTION,  PHASE  4 

1  INTRODUCTION 

This  report  documents  the  results  of  Phase  4  of  the  Terrain  Data  Extraction  study. 
Phase  4  effort  was  performed  during  the  period  September  1982  to  November  1983. 

The  overall  study  addresses  the  problem  of  extraction,  from  digital  aerial  imagery, 
of  terrain  features  essential  for  decisions  associated  with  military  operations 
such  as  cross-country  movement.  Also  addressed  in  the  study  are  some  aspects  of 
the  generation  of  terrain  analysis  products  (e.g.,  maps)  necessary  for  the  deci¬ 
sion-making  associated  with  these  military  operations.  Primary  effort  in  the 
overall  study  prior  to  Phase  4  has  focused  on  the  extraction,  from  panchromatic 
aerial  photography,  of  vegetation  and  other  land-cover  boundaries,  and  on  the 
extraction  of  forest  features  such  as  canopy  closure  and  stem  spacing.  See  the 
Phase  1  report  ETL-0241 ,  November  1980  for  a  more  complete  description. 

Phase  4  involved  four  tasks: 

(1)  Produce  automated  delineation  (boundary  outlines),  based  on  extraction 

techniques  developed  in  Phases  1-3,  of: 

o  Vegetation  classes 
o  Percent  canopy  closure  intervals 
o  Tree  stem  spacing  intervals 

(2)  Develop  algorithms  and  provide  software  for  scaling  and  mosaicking 
vegetation  data  derived  from  separate  adjacent  scenes,  as  appropriate 
for  covering  a  given  map  area. 

(3)  Develop  algorithms  and  provide  software  for  converting  area  boundary 
(outline)  data  to  a  format  suitable  for  input  to  the  ETL  Interactive 
Graphic  Design  System  for  use  in  polygon  manipulation  programs. 


(4)  Provide  a  detailed  design  for  an  interactive  system  capable  of  performing 
the  tasks  (1),  (2)  and  (3),  above. 

Task  (4)  was  completed  early  in  Phase  4  and  documented  in  a  report  submitted  to 
USAETL  dated  January  1983  and  entitled  Interactive  Image  Analysis  System  Design. 

The  results  of  Tasks  (1),  (2)  and  (3)  are  documented  in  this  report  in  Sections 
2,  3  and  4  which  follow.  Software  developed  in  these  tasks  is  not  included  in 
this  report  but  has  been  delivered  separately  to  USAETL.  In  Task  (3)  not  only  was 
the  polygon  conversion  problem  addressed,  but  attention  was  given  also  to  the  auto¬ 
matic  extraction  of  polygon  boundaries  for  a  complex  theme  representing,  for  example 
all  forest  areas  in  a  large  scene  (i.e.,  multiple  polygons,  polygons  with  holes, 
polygons  within  polygons,  etc.). 


DELINEATION  OF  VEGETATION  CLASSES,  CANOPY  CLOSURE  AND  STEM  SPACING 


Image  products  depicting  automated  delineation  of  the  following  have  been 
produced  using  extraction  techniques  which  have  been  developed  in  this  study: 

•  Boundary  Outlines  of  Vegetation  Classes 

•  Boundary  Outline  of  Forest  Canopy  Theme  According  to 
Percent  Canopy  Closure 

•  Boundary  Outline  of  Forest  Areas  According  to  Average 
Tree  Stem  Spacing. 

Figure 2-1  shows  the  basic  digital  aerial  image  from  which  the  above  vegetation 
feature  outlines  were  extracted. 


Figure  2-1.  Digital  Aerial  Image 


Vegetation  themes  were  produced  from  the  basic  image  using  techniques  developed 
in  this  project  and  described  previously.  In  this  particular  scene  only  three 
broad  classes  exist:  Forest,  Grass  and  Bare  (includes  roads  and  buildings). 
These  themes  have  been  color  coded  in  blue,  cyan  and  white  respectively,  and 
superimposed  on  the  basic  image.  See  Figure  2-2. 


The  result  is  a  gray-level  canopy  closure  map  in  which  pixel  brightness  is  pro¬ 
portional  to  percent  canopy  closure.  In  turn,  that  map  is  level  sliced  into 
five  equal  pixel  brightness  intervals,  and  a  color-coded  theme  is  produced  for 
each  interval.  Figure  2-3  shows  these  themes  superimposed  on  the  original  digital 
image.  Color  coding  is  as  follows: 


Canopy  Closure  Interval  Color  Code 


0  -  20% 

20  -  40% 
40  -  60% 
60  -  80% 
80  -  100% 


Red  (Reddish-gray  when  super 
imposed  on  image) 

Cyan 

Blue 

Green 

Purple 


Figure  2-3.  Scene  With  Canopy  Closure  Themes  Superimposed 


Tree  stem  density  themes  are  produced  from  the  basic  digital  image  by  first 
producing  a  binary  map  in  which  single  pixels  representing  tree  crown  locations 
are  given  a  value  of  "1",  with  "0"  for  the  remainder  of  the  image.  This  crown 
location  theme  is  extracted  from  the  basic  image  using  the  crown  template  and 
processing  procedure  developed  and  described  previously  in  this  study.  Next, 
the  crown  (stem)  theme  is  modified  such  that  each  crown  location  pixel  becomes 
a  3  x  3  pixel  array,  with  each  pixel  having  an  intensity  value  of  163.  This 
modified  theme  is  now  subjected  to  a  27  x  27  pixel  (729  pixels)  moving  average. 
A  gray  level  image  results  (see  Figure  2-4)  in  which  pixel  brightness  is  propor¬ 
tional  to  stem  or  crown  density.  Since  27  x  27  pixels  correspond  to  22.1  x 
22.1  meters2  in  this  image,  a  pixel  brightness  of  |^x'Xf7Six~e'ls  ^or  ^ 
corresponds  to  an  area  per  tree  of  22.1  x  22.1  =  488  meters2  . 


Figure  2-4.  Moving  Average  Applied  to  Tree  Crown  Theme 


The  gray  level  image  resulting  from  the  27  x  27  pixel  moving  average  is  level 
sliced  at  appropriate  digital  levels  so  as  to  produce  themes  representing  in¬ 
tervals  of  average  tree  area  or  stem  spacing.  These  themes,  plus  the  3x3 
pixel  crown  locations,  are  shown  in  Figure  2-5. Theme  color  codes,  and  the  stem 


Figure  2-5.  Average  Stem  Spacing  Themes  With  Crown  Location  Theme 
spacing  intervals  which  they  represent,  are  as  follows: 


Theme  Color  Code 


Stem  Spacing  Interval 


yvv  V  V  VVV  VVVV  V  v  .• 


Figure  2-6.  Scene  With  Stem  Spacing  Interval  Themes  Superimposed 

spacing  vs  digital  level  number  in  the  gray  level  image  of  the  averaged  crown 

location  theme,  it  is  assumed  that  average  tree  area  is  inversely  proportional 

k  Q77 

to  digital  number,  or  A  =  Since  DN  =  2  corresponds  to  488M2  ,  A  =  . 

Using  a  hexagonal  tree  model  (see  Figure  2-7)  it  can  be  shown  that  stem  spacing 
S  =  2h  =  1.075A^.  Thus  S  =  33.6  (DN)"^.  Average  stem  spacing  for  five  selected 
digital  level  numbers  are  tabulated  in  Figure  2-7.  The  figure  also  shows  the  loca¬ 
tions  of  those  digital  levels  in  the  histogram  of  the  image  of  the  averaged  crown 
location  theme. 

Finally,  the  three  digital  images  of  themes  for  vegetation  classes,  canopy  clo¬ 
sure  intervals  and  stem  spacing  intervals  were  processed  to  produce  new  themes 
representing  the  boundaries  of  the  original  polygon  themes.  These  polygon  boun¬ 
dary  themes  have  been  converted  to  the  SIF  format  and  loaded  onto  three  digital 


3  IMAGE  MOSAICKING 
3.1  INTRODUCTION 

A  mosaic  is  an  assembly  of  individual  photographs  fitted  togetner 
systematically  to  form  a  composite  view  of  the  entire  area  covered 
by  the  photographs.  The  mosaic  gives  the  appearance  of  a  single 
photograph.  producing  a  complete  record  of  the  area  photographed 
Perfect  matching  of  the  images  of  adjoining  photographs  is 
virtually  impossible  because  of  variations  in  altitude  of  the 
airplane,  large  differences  in  elevation  of  the  ground,  divergence 
of  the  camera  ails  from  the  vertical  at  tr.e  time  of  exposure,  and 
errors  inherent  in  the  photographic  system  of  camera.  lens.  film, 
print,  paper  , etc.  By  the  use  of  proper  techniques,  however,  these 
errors  can  be  held  to  a  minimum  and  not  permitted  to  accumulate 

Tc  construct  an  accurate  mosaic,  all  individual  photographs  must 
be  reduced  to  a  common  scale,  must  be  rectified  to  remove  tne 
distorting  effects  of  camera  tilt.  and  must  be  assembled  tc 
accurate  ground  control  at  selected  photo-identifiable  points  for 
control  of  azimuth,  scale,  and  position. 


Mosaics  are  very  useful  in  all  types  of  planning  activities  and 
have  the  advantage  of  showing  as  a  single  picture  a  wealth  of 
detail  of  the  entire  area  under  study. 


3.2  DISPLACEMENT  OF  IMAGES  DUE  TO  TILT 


At  a  result  of  tilt  alone>  the  images  appear  to  be  displaced 
radially  toward  the  isocenter  on  the  upper  side  of  the  photograph 
and  radially  outward  or  away  from  the  isocenter  on  the  lower  side. 
Along  the  isometric  parallel  (line  through  the  isocenter 
perpendicular  to  the  direction  of  tilt)  there  is  no  displacement 
relative  to  an  equivalent  untilted  photograph 


Ur>4i  Phoio  •  Plane  r 


Figure  3-1.  Image  displacement  due  to  tilt 


The  Figure  3-1  illustrates  the  geometry  of  an  aerial  photograph 
and  the  symbols  have  the  following  interpretation: 

O  *  Location  of  camera  lens 

OP  =  Optical  axis  which  is  perpendicular  to  the  plane  of  the 
photograph  at  the  principal  point  P  (geometric  center > 
This  is  equal  to  the  focal  length  of  the  camera. 

t  =Tilt  angle 

i  =  Isocenter 

n  “Nadir 

Di  “  Distance  of  image  from  the  isocenter  along  the  principal 
line 

f  =  Camera  focal  length 

dt  “  Displacement  of  image  relative  to  the  equivalent 
untilted  photograph 

Then,  the  magnitude  of  displacement  dt  is  given  bu  the 
equation 

3. 

DL 


dt  =  - - — 

(  -f/s.nCi)  -  Di) 


when  t  is  very  small  and  expressed  in  radians. 


3.3  PROJECTIVE  TRANSFORMATION  EQUATIONS 


Assuming  that  light  rags  travel  in  straight  lines*  that  all  the 
rags  entering  a  camera  lens  sgstem  pass  through  a  single  point  and 
that  the  lens  sgstem  is  distortionless*  then  a  projective 
relationship  exists  between  the  photographic  coordinates  of  the 
image  points  and  the  ground  coordinates  of  the  corresponding 
object  points  as  illustrated  in  Figure  3-2. 


r.  \h 
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Figure  3-2.  Exterior  Orientation 


where 


X,  Y,  Z 

Xj<  Yj,Zj 

v  c  c 
Xi,  Y i .  Zi 

x,  y,  z 

w,  4-  k 


Three-dimensional  rectangular  coordinate  system  *cr 
the  object  spacr 

The  coordinates  of  any  point  J  in  the  object  space 
The  position  of  the  exposure  center  of  photograph  1 
The  photo- coor d i nate  system. 

Rotation  angles  about  7»y,T  respectively  which 
define  the  direction  of  the  optical  axis.  The 
angles  are  positive  in  counter-clockwise  direction 
When  w=<|=k=0.  then  x»y»7  frame  is  parallel  to  the 
X, Y, Z  frame  and  the  optical  axis  is  perpendicular 
to  the  X-Y  plane. 


Then*  the  projective  tT ansf or m^r ion  equations  can  be  written  a  own 

as 


ux 

I 

•  — l 

X 

*  ^'j  j 

.  ™,r  Cy.j  -»  ^  m„  C-f)] 

Yj  -  Y* 

5  *ii  [ 

^ x Jj * *t  ma*  Oj-yp)  *«• 

Z-j  - 

* 

c*ij-xp)^  ^xs  (yx-xp)T  m3aC-‘pi] 

where 


^ii  r  c^sC^ijCosCk;) 

TnJX  *  cos  Cwt)  .Str*  C  k‘.)  f  s'ir>  ( u>i)s in Oc'*) 
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">15  *  «m£o;)  StnC*')  -  Ces(Yo.).  sin  ($))  •  Cas  (fc.) 

n*.«  a.  -  c>&*  C&\)-s\n(K>) 

maa  *  Cos£u\).  CosCfc.)  -  fcmCwO-siflC^.)  «*'n(V.) 

*  Sin  •  <*>s(k.)  —  CosCw't)  •  sfn  £♦*,)  *  s,r>  C^)i 

IO31  s  sin  (_^i) 

rn 33.  *  -  sin  £co.)  ■  c*sC4!) 

10^-5  «  Co<  C  to.)  •  O sC*0 

*jj»  y;j  *  Image  coordinates  of  an  object  point  J  on  photo  i 
>ij  *  The  photo  scale  factor  at  image  point  J  on  photo  i 
=  in  the  newly  rotated  photo-coordinate  system  about 
=  the  T.y ,~i  axis  of  the  photo. 

When  w*t~k*0>  then  referring  to  photo-coordinate  system  as  Fe .  y"0 
the  following  relationship  exists: 

Xj  -  X?  *  >ij  ** 

Xi  ”  Y;  =■  >J  5. 

c 

£« 

In  this  case.  x0  .  ye  <  and  zc  are  parallel  to  X.  V.  and 


respectively. 


the  lower  pert  of  the  principal  line  which  passes  through  the 
nadir  point  (n).  The  direction  of  tilt  with  respect  to  the  ground 
reference  coordinate  system  <X,Y»Z>  is  defined  by  the  azimuth 
angle  <«<)  which  is  measured  in  the  X-Y  plane  clockwise  from  the 


positive  Y  axis 

to  the  projection 

Of 

the  X-Y  plane. 

The  elements 

in 

the  projective  transformation  equation 

expressed  in 

terms  of  t,  s.  and^as 

follows: 

mu 

9 

-  c.**C*3‘  c»sC«0 

— 

*>tn  C«-)  tin  (•*) 

■x. 

Si  V)  C*  )  c,*a 

- 

Co  S  Cs)  •  Cos  (.■<•)  si"  («0 

m 

—  Sin  (•*)  Sin  (.<) 

nrta, 

■e 

C«S  (s)  ■  SinO) 

- 

sin  (s)  •  C-stt)  .  CoS  U) 

•"U 

= 

-  SinCs3-  sin£»0 

- 

CjosQs)  CosC*-)  -  CO  sO) 

• 

-  C os  C«<)  •  S«*n  Ct) 

) 

-  Sin  Cs)  Sin  Ct) 

- 

*  CoS  Cs}  •  Sin  Ct ) 

r 

c#s  tt; 

A  major  disadvantage  of  this  system  of  defining  orientation  is 
that  it  breaks  down  for  the  ideal  aerial  photograph  which  is  truly 
vertical.  When  there  is  no  tilt*  the  two  planes  are  either 
parallel  or  coincident*  and  there  is  no  line  of  intersection  (or 
principal  line);  thus*  the  angles  of  swing  and  azimuth  are 


undefined. 


3.3.2TRANSFORMATION  MATRIX 


The  projective  transformation  equations  can  be  written  in  matrix 
form  as  follows: 


1 

z 

» 

—  o 

_ 1 

T 

*‘j  -xr 

T 

**»» 

«r>,a 

vj  -xc 

-  "V » 

*ii  - 

mxa. 

mis 

c 

zi  -Zi 

-f 

mu 

tTji 

_  _ 

, 

Then. 


—  — 

—  — 

®ii  *n |2  (itj 

| 

••'ll  Hii  i»ii 

fa  Vj  Zjj  =  [xf  rf  z?]--x,j  l>r  vrfj 

Flat  s»i  (ijj 

jrj  0] 

mu  ('ll  •Hij 

(03*,  rjj 

m3i  m»a  e»s? 

1 

— 

or 


o  n 


xHX 


Where 

«  a  vector  of  points  in  the  object  space  - 
«  a  translation  vector  for  the  principal  point 
*  a  vector  for  image  points  in  the  photo-coordinate 
system 

C  R  3  *  a  rotation  matrix 

'Xjj  *  an  overall  scale  factor  due  to  camera  altitude 
Using  the  homogeneous  coordinate  system*  this  can  be  written  as 

“  v  z  n  Ajc,  v  -f  n 

3.33  ESTIMATING  THE  TRANSFORMATION  MATRIX  CTRFI 

From  the  transf ormation  matrix  derived  above*  the  equations 
relating  an  object  space  point  to  an  image  point  are 

X  =  *  4-  ma,  }  +  ) 

y  =  (  x  +  max  y  4* 

Z  s  ^  Vf  \ | •,  x  4-  (n  ^  4  ^  Tx  ” 
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The  overall  scale  factorVij  is  altitude  dependent  and  is  generally 
known.  Then.  the  remaining  nine  unknowns  can  be  determined  by 
using  three  control  points.  If.  however,  more  than  three  control 
points  are  available,  then  it  becomes  an  over  determined  system  and 
a  least  squares  technique  may  have  to  be  used  The  least  square 
problem  is  formulated  as  follows:  These  equations  have  a 

general  form 

UJ  =  a  tV4+  cv2 

where 

W  is  a  dependent  variable 

Vj  ,  V2are  ind  ep  t-;;d  ent  variables 

a.  b.  and  c  are  the  parameters  to  be  determined 

Then  the  normal  equations  to  be  solved  for  the  least  squares 
estimates  of  a.  b.  and  c  are 


Zw  * 

no. 

-Y 

c  Z  vz 

1  WVi  - 

aZ  Vt 

■T 

c  Z  vt  v 

IWVs  * 

Q  2.vz 

-t 

b  2  V,v2  4 

C  2.  V* 

where  all  the  summations  are  over  the  number  of  observations  (n). 

Thus,  the  resulting  nine  equations  when  solved  provide  estimates 
for  the  nine  elements  of  the  transf ormat ion  matrix  CTRFD.  Using 
this  matrix  CTRFD.  a  point  operation  can  be  performed  for  each 


image  point  to  generate  its  object  space  equivalence.  However, 
this  way  of  doing  rectification  is  very  time  consuming  and 
expensive. 

Shortly,  an  analytical  method  will  be  outlined  which  consists  of  a 
conceptually  simple  model  of  the  exterior  orientation  of  an  aerial 
photograph.  The  necessity  of  such  a  model  is  evident  from  the 
fact  that  a  set  of  fast  and  efficient  operations  are  needed  to 
carry  out  image  rectification. 

The  algorithm  that  performs  a  least  squares  fit  to  the  control 
points  is  implemented  in  the  program  MOSAIC. 


•bed  is  the  user  specified  output  window 

ABCD  is  the  input  image  window  to  be  extracted  * 


Now.  the  transformation  matrix  CTRF3  maps  ima$'e  points  to  the 
corresponding  object  space  points  as  follows: 

[X  Y  2]  »  "Vij  C  x  y  13CTRF3 

Then. 

Cx  y  13  *  •  CX  Y  ZHTRF3  * 

This  enables  computation  of  image  coordinates  corresponding  to  the 
user  specified  output  window  coordinates. 


Now.  if  the  user  specified  output  window  "abed"  has  the  map 
coordinates  (Xa.Ya).  (Xb.Yb).  (Xc.Yc).  and  (Xd.Yd),  then  the 
corresponding  image  coordinates  given  by  the  above  equation  are  <  xa 
.yA).  <xk.  yb>«  <*c*  yc>'  *n<*  <xa'  Me*  respectively.  However,  the 
input  image  window  "ABCD”  to  be  extracted  for  mosaicking  is  given 


MIN  <xA, 
MIN  <yA. 


For  upper  left  corner  of  window 


and 


*  MAX  <  xA ,  xfe.  xe, 
yc  -  MAX  <ya.  yfc,  yc, 


For  lower  right  corner  of  window 


This  algorithm  is  implemented  in  the  routine  IOWNDO. 
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V  W  -p  v 


p 

I 

N 


%,  y 


A 

* 

4>o 


«*  «r  r.  -  -  ,  -r 


•  ;  •  .  <r,  •  *.  • .  * 


■  Principal  point 

*  Itocenter 
«  Nadir 

■  Orthographic  projection  of  X, Y  on  a  plane 
parallel  to  X-V  plane  with  origin  at  P 

*  Translation  vector  with  components  Ca^ayl 

■  Angle  between  the  photo  coordinate  system 
Cxr.yrD  and  the  rectangular  system  <x.y> 

*  Angle  between  the  pricipal  line  and  xp. 


Then. 


X 

-1 

—  — 

-  a. 

•r  s 

R  Co) 

Y 

1 

zs 

< 

1 

A. 

>* 

_ 1 

where 

R  1((j))  *=  Rotation  of  image  data  from  (x  .y  )  system 

to  <  x. y )  system 

M  *  Scale  matrix  which  removes  the  tilt  effect 

s  *  Overall  scale  factor  due  to  altitude 


kx.Ay  ■  Coordinates  of  the  isocenter  with  respect 

to  the  principal  point 

It  should  be  recognized  that  <&x.&y>  are  dependent  upon  the  tilt 
and  swing  angles  of  the  photograph  which  are  not  known  to  begin 
with.  Hence.  as  explained  later.  they  are  assumed  to  be 
negligible  initially  and  this  model  enables  their  determination  by 
using  photo-coordinates.  Later,  an  iterative  improvement  is  made 


6 


to  the  predicted  values  of  tilt  and  suing 


Then#  writing  in  terms  of  x ;  .  y.  (i.e.,  image  coordinates  wit 

respect  to  isocenter  <1)  as  origin)#  the  model  is  stated  as 


Since#  LMj  is  a  diagonal  matrix.  ar.c  [Ml  are  commutative 

that  is 

Cr'cwIL^i  -  tMlt *"<»] 


therefore,  one  can  wnte 


R  it  a  point  in  the  object  space  with  spherical 

coordinates  <r,$,  •) 

f  is  the  focal  length  of  camera 

t  is  tilt  angle  and  it  is  around  x'axis 

(i.e.,  in  the  Y'-Z'plane) 

R*  is  the  projection  of  R  on  the  horizontal  plane 

(i.e.,  parallel  to  x'-y' plane)  at  distance 
f  from  the  origin.  The  spherical  coordinates 
of  R1  are  denoted  as  (r'»4'9>- 
R°  is  the  projection  of  R  on  the  tilted  photo  plane. 

Its  spherical  coordinates  are  denoted  as  <r",^,  ft). 
P.l.N  are  the  principal  point,  isocenter,  and  nadir 

respectivelg. 

Then,  the  vector  R'  can  be  written  down  as 


and 


But 


^  ■  r  C  «  •*  j  •m(ejcin(l)  ♦  kCci(,()] 


o"  l.  /  -  .  — 

K  •  *  V.  *  ««"(.«)  C«s(.4»)  j  S'"  (W  «'•<(«)  +kc<»sC»)) 


**’  •  *  -  f  -  r'  C.*(p) 


(3.3) 


•  •  • 

Since,  the  tilt  is  around  X  axis,  OP  is  in  the  Y  -Z  plane  and  the 
principal  line  is  along  Y*axis.  A  unit  vector  in  the  direction  of 
OP  is 


f>  •  j  t'm  (t)  4  C  C*<  C*-) 
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f?  .  n 


OP  -  f 


r"  C  *»"03  •  ft*"  CP)  •  6i*(.U  * 


Al! 


Therefore 


equating 


(3.3)  to  (3.4) 


Y  *  C~OS  (0) 

»  y" 

Cst'nC©)-  Sin (£>)•  sin  (t)  ■+• 

cos  (jo)  -  cos(t;  } 

• 

r 

a 

=  r 

(  s in  £fl)  .  sinfcfi)  *  -L 
V  Co*  CO) 

-4- 

c-o e  C6>  .  fi-  if  ^ 
cos  C©)  ^ 

*  if 

v  —  v 

*We)  .  -k.nf^.l 

t* 

V” 

Co sC©) 

z 

= 

v  */*.• 

j.  ' 
~  I 

)•*• 

=  ( 

Jl  S‘*n^)-6  in  (dp)  — 
<*f 

GI 

~r 

)<■ 

£ 

C  -  ^0-t 

* 

IRi  4. 

A  r 

T7 

S 

rp;  . 

-r  • b 

At 

r' 

» 

w  # 

JL  .  15.  t 

Y'  f 

But  angle 

»  n 

IRe  I  - 

t  and  angle  I^IRo® 

<  TT 

-  t>/2 

therefore. 

angle 

u"r>  err  -  t  >/2 

and  IRo=  l"Ro 

Hence. 


At 

“ 

which  implies 


Dl  .  t 

“7 
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•-v'  aJi aC aSkjjJCi 


Ax 

X# 


A  y 

Y 


Define 


Then, 


<Mx  -  1) 


—  and  (My  -  1) 


(Mx  -  1)  *  (My  -  1 )  =  Di.  t/f  or 
Mx  *  My  =  (1  +  Di. t/f) 


3k 5.2  PREDICTION  OF  PHOTO  PARAMETERS 


Referring  to  Figure  3-5,  for  a  point  (xr,yr),  the  distance  from 
the  isocenter  along  the  principal  line  is  given  by 


Di  ■  Xp.  cos(^o)  +  y p  sin(l^o)  +  f.  t/2  and 


My*  My  *  1  *+• 


(  xi> Cp.}  *  y t. 


j  .  41  .  y^sinCA,))t 

~  ***  T  « 


Substituting  for  Mx  and  My 


— -C->sC^^^*ln  Cd>>^  ^  |(*i  -  y»  (3.5) 

CX-Qy)*  S  C  tejjt  >_p  t  j.  (» ;  5~.  *(»)  4.  y ;  Cq«(»))  (3.6) 


The  unknowns  in  the  above  equations  are 


**  »  a  y 


,  t  -  <t>o  >  as  the  overall  scale  factor  s 


and  the  focal  length  f  are  generally  known. 


3.5.2. 1  TRANSLATION  VECTOR  7T 


The  two  equations  in  above  provide  the  following  relationship 


*  !  C.O*-  —  y  \ii<r\£s>) 

X*  *in  C4*l  y i  C,o&(d>) 


*•  y  i  -4  dr)  (a> J 
X;  y; 


Sc-i vino  for 

_  (■Y-°x)xj  —  Cx-Q,)y; 

.+.  (Y^dy)yi 

which  involves  ti.ree  ur<  i>  nouin  • 

I  3  y  I 


Using  three  control  points.  ;ai.  C^l.  CP2.  and  CP3i  these  unknown* 
can  be  determined  as  follows 


Denote  by  <Xj.»Yi).  <Xa.>Yi)i  and  (XS>Y5)  the  map  coordinates  of  the 
t  h  t  e  e  control  points  and  to  (  xA  »  y1  ) ,  **»'¥***  and  <x,<y*>  the 
corresponding  image  points  respectively.  Then,  using  two  control 
points  CPI  ,  and  CPS 


•W(*>  -  -  C 

(  Xi-fl*)*!.'  (.X"  fl0  *4' 

c  r  e  4  »  “tn  u  1 1  i  p  lying. 


C^a.  ~Qy)  **>  ~  Cx*~a  ►■)>*.' 
C  Xi -a,)**;  ♦  ( 


LC  x-« 


*±t  “C*!  -fl*)  JfiQ  £C  *>■"**)  C  ta"  >afj 


ol  point1 


Cy  t>an*4'<  ng 

*■1 

a  S  y  ^  X 

n,  (*»-*.) 

“  l  c-v^‘ 

C'nr-'n.) 

•{'"U&l-'-l 

o-* 

* 

1 

o-> 

«i 

*V*n* '  \ 

u 


\VVWV  "WW  \ 


=  R 


Then  or.r  get-. 


P  A*  •+  ^  Q*  ■+  R.  «:  O 


This  quadratic,  in  a  „  has  twc  roots 


^  -q  ±  yo‘-^pR 

C4,i  *  - - - - - 

IP 

Since.  a%  and  ay  art-  components  of  a  trar,aatior.  vector  A.  they 
should  be  real  numbers.  Also.  if  the  origin  of  the  map 
coordinates  system  is  selected  such  that  the  photograph  falls  in 
the  first  quadrant,  then  both  ak  and  ay  must  be  positive  and  a 
correct  pair  can  be  always  be  selected. 


3. 5. 2. 2  ROTATION  ANGLE  {> 


One*  a,  and  ay  ara  known,  the  angle  can  easily  be  determined  by 
evaluating  the  following  equation  for  control  point  CPI: 

-I  0*(4?)  n  *i.1  -  C*i-*0  y kl 

(X.-a.)**.-  -+  CYi-a.)  >*; 

Then#  the  rotation  angle  is  determined  such  that 

A 

-  HL  <  $  <.  it 

3 .  —  t  —  a, 

3. 5.2.3  DIRECTION  OF  TILT  <$o> 

Equation  3.1  is  re-written  as 


X  =  H-  S  M*  •  £  *;  cosC«p)  - 


Substituting  for  Mx 


L  (  xpCog(f,)  4-  yi»sin(<gj 

1  l  * 


T 


X  -  «> 


o  ;  cos(&)  -  am 


Using  control  point  CPI 


.  ^  i*  .  (  «t  -»  >jSin«.,)  •) 

1  C  f  j 


Xi  "  Q. 


s  C*i«  Os»)-v4;aink») 


-P.  0.7) 


Similary#  using  control  points  CP2  and  CP3 


**  *  t 


*a  CosCffaW  3*  s»n(q>.  )  j  ^ 


i 


Xj  • 


»  (x»',  cos  t*J  -  Ja;  sin  »>)) 


;P*  (3.8) 


3. 5. 2.4  MAGNITUDE  OF  TILT  (t) 


Now,  equation  (3.10)  is  solved  for  tilt  <t>  as  follows: 

t  .  f  c  t-i  -  ft.) _ 

C*« -*  s)  «■  u.- y*) 

This  provides  tilt  (t)  as  a  signed  number.  This  signed  value  of 
tilt  (t)  with  the  signed  value  of  (^o  permits  computation  of  swing 
angle. 

If  t  >  0  ,  then 

Swing  angle  -  3TT/2  -  (jk>  and. 

If  t  <  O  i  then 

Swing  angle  =  ^ /2  -  $o  . 


This  algorithm  is  implemented  in  the  routine  RECTFY. 


3.5.3  ITERATIVE  IMPROVEMENT  OF  PHOTO  PARAMETERS 

The  first  estimation  of  photo  parameters  is  made  by  assuming  the 
separation  between  the  principal  point  and  isocenter  as  negligible 
because  it  depends  upon  the  photo  parameters  to  be  predicted.  The 
first  prediction  of  the  tilt  and  swing  angles  under  that 
assumption  allows  computation  of  the  isocenter  coordinates  (a* ,ay) 
which  can  be  used  for  the  next  iteration.  The  revised 
prediction  of  the  photo  parameters  resulting  from  such  iterations 
shows  a  significant  improvement.  Re-stating  the  basic  model 


*  -  An 

s  s  •  0t>>  W 

't  af  Yr  - 


U*ing  the  vector  notation  it  can  written  a: 


x  -  A  +  5[r'w)][m][^-i] 

A  r  X  -  s  [RM)][m]v  +  s[RW] 

=  x  -  (►-«) 

Ui  h  i-  r  e 


^  =  s.[RC*)][M]4 

Now.  the  first  predictions  are 

A*  ,  £  ,  *.  ,  t 

and  they  enable  computation  of  rn*3.  and  such  that 

\  «  *  L  p'(V)][m#]  *r 


satisfying 


Then<  the  iterations  can  be  stated  as 


Ax  -  x-a-^o  * 

neto  vm-loas  op  >  i>  >  [M  ]  A 

^  =  s  •[  r'^^TQm1]  A1 

Ai  -  X-C?t-5t)  -  A*  -Cfi-'fe-S,) 

Cc-V  n«J  M»|«ts  eA-  <£  ,  ^  ^  )  t  ^  ^ 

0n<^  re^aah  4he  i Ve*r«Ucm  . 

and  the  convergence  to  be  tested  is 

(  A*  -  A**'  j  <  i  t>',xd 

Thus,  once  the  position  of  the  isocenter  is  determined  within  a 
subpiiel  accuracy.  the  iterations  are  stopped.  The  resulting 
values  of  the  photo  parameters  are  used  in  the  subsequent  image 
rectification  procedure.  During  algorithm  testing,  it  has  been 
observed  that  the  convergence  is  reached  within,  generally,  a  couple 
of  iterations. 
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S*.S 
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3.6  IMAGE  RECTIFICATION 


Assuming  that  the  distortion  doe  to  film  shrinkage  during  film 
processing  is  negligible,  a r.  aerial  photograph  can  be  considered 
essentially,  rectified  when  the  distortion  due  to  tilt  effect  is 
removed 


It  is  evident  that  the  scale  matrix  CM3  removes  the  distortion  due 
to  ! lit  Nou. 

O 

O  My 

0»<*  M*  -  My  =  1  «(*»)+» t 


wh  ere 


t 

f 

4>o 


^  *  f 1  yp  ^ 


tilt 

camera  focal  length 

rotation  angle  between  the  principal  line 
and  the  x  direction  of  the  photo-coord  inatc 
sy s  y tm. 

photo -coord  mates  of  an  image  point 


The  quantity  (  Xp.  cos  )  +  yp  sin((^o)  )  is  a  distance  along  the 
principal  line  of  an  image  point  <xp  . y^ ) .  Since,  the  pnoto 
digitization  may  occur  in  some  arbitrary  direction.  it  mag  be 
inefficient  to  apply  the  scale  matrix  CM]  to  image  data  keeping 
the  scan  direction  unchanged  This  is  due  to  the  fact  that  fo" 


each  pixel  on  a  scan  line  the  distance  along  the  principal  line 
mag  be  different/  implying/  that  it  will  have  to  be'  computed  for 
each  pixel  on  a  line  and  for  each  line. 

It  is  the  intent  of  this  section  to  show  that  there  is  a  preferred 
scan  direction  for  each  photograph  so  that  the  scale  matrix  Chi 
can  be  applied  to  remove  the  tilt  induced  distortion  as  a  fast 
scan  line  oriented  operation.  This  is  illustrated  in  Figure  3-7. 


Vp 


Figure  3-7.  Actual  and  preferred  scan  directions 


for  an  aerial  photograph  digitization 


In  the  Figure  3-7  the  rectarjoiar  p h ot o-c cor d inat e  system  <Xp.vp> 


has  an  origin  at  F,  the  principal  point.  which  is  also  the 
geometric  center  of  the  photoglyph;  the  o  •  incipal  line  Pr  and 
isometric  parallel  Is  form  a  rectangular  coordinate  system  with 
origin  at  I.  the  lsocenter;  and  the  two  coordinate  systems  are  at 
an  angle  <^o.  The  photograph  u  digitized  along  direction  "Ail  so 
that  the  scan  lines  are  p a r  a  1 1 e 1  to  Xr  -a*is  The  pre+errec 
direction  of  digitization,  in  this  case,  i:  "ab"  which  is  parallel 
to  the  isometric  parallel  line  Then,  it  is  easy  to  see  that  for 
each  scan  line,  all  pixels  or  the  line  will  have  the  same  distance 
along  the  pncipai  lir. 
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where  Li  is  distance  of  an  imaos  point  from  the  isocenter 

I  along  the  principal  line  Pr 

Suppose  that  the  image  data  :  =  rotated  to  achieve  the  preferred 
scan  direction  Tr.en.  Figure  3-8  illustrates  the  rectification 


process  as  a  line  oriented  operation 


x-  component— 


Figure  3-8.  Image  rectification  concept 

Let  (x^.gr>  be  the  coordinates  of  an  image  point  in  Xp-Yp  frame 
with  R  as  the  origin,-  <*,-.«;)  be  the  coordinates  of  the 

•ame  image  point  when  the  origin  is  shifted  from  P  to  I. 

Clearly, 


Now,  for  the  image  points  J,  K,  and  L  on  the  preferred  scan 


direction  "ab",  their  corrected  positions  are  J',K'  and  L' 
respectively  which  are  radially  outward  from  the  isocenter  I. 

When  the  vectors  JJ'  ,  KK '  »  and  LL '  are  represented  as  two 
orthogonal  components,  one  perpendicuar  to  Pr  and  the  other 
parallel  to  Pr,  the  latter  components  turn  out  to  be  the  same  for 
all  points.  The  component  perpendicular  to  Pr  is  the  x-COmponent  in 
the  rotated  system  and  by  the  same  token,  the  component  parallel 
to  Pr  is  the  Y-Component.  Also,  it  is  known  that 

X-  «  (Mk)  X,-  =  +  *. 


Y  -  Um  e  (M^)  y; 

are  coordinates  of  the  corrected  positions  of  image  points  in  the 
x;  _y(  system.  But  when  the  image  is  rotated  to  get  "ab"  as  the 
scan  direction,  it  is  in  the  Is-Pr  coordinate  system  such  that 
each  scan  line  is  a  constant  distance  a way  from  the  isocenter  I, 
implying  the  Y-component  for  all  pixels  on  line  "ab"  is  the  same. 
This  means,  as  a  result  of  rectification,  the  line  "ab"  has  been 
shifted  from  its  previous  location  in  the  rotated  system.  The 
X-components,  at  the  same  time,  are  within  scan  shifts  as 
illustrated  in  Figure  3-8.  Thus,  the  rectification  is 
accomplished  as  a  line  oriented  operation  in  terms  of  a  constant 
shift  of  the  input  line  and  a  position  dependent  shift  for  each 
pixel  of  the  input  line.  Also,  it  should  be  noted  that  a  position 
dependent  shift  for  pixels  is  proportional  for  all  lines,  meaning, 


a  table  can  be  set  up  only  once  and  applied  to  all  lines.  The 
advantage  of  rotating  image  data  to  achieve  the  preferred  scan 
direction  is  obvious  now.  It  enables  implementation  of  fast  line 
oriented  operations  to  accomplish  rectification  due  to  tilt 
induced  distortions.  Symbolically,  if  ®<  and  are  image  data 
coordinates  in  the  rotated  Is-Pr  system,  then  after 
rectif ication,  the  new  image  data  coordinates  will  be  e^g  and  (3  c 
such  that 


o<c  a 

« 


C  i  +  po  < 
Cl  -  p)  p 


where  sign  is  that  of  the  angle  Qo. 


When  image  data  are  rotated  to  achieve  the  preferred  scan 
direction<  the  isocenter  (I)  becomes  the  origin  of  the  rotated 
coordinate  system  <Is-Pr>.  The  rectification  takes  place  in  this 
rotated  coordinate  system.  The  coordinates  of  (I)  in  terms  of 
lines  and  pixels  is  derived  as  follows: 

If  the  input  image  window  is  m  x  n  (i.e.,  m  pixels  and 
n  lines),  then  the  isocenter  coordinates  with  respect  to 
a  are  given  by 

Ipixel  =  0  5<m.  sin<  (|>o  )  +  n.  cost  (^o  )  ) 

Iline  =  0.  5(m.  cost  )  +  n.  sin(  0o  )  +  f.  t  ) 

These  coordinates  permit  computation  of  distance  Di  for  each  scan 
line  and  pixel  in  the  rotated  system. 


3.6.2  FACTORING  ROTATION  INTO  TWO  SEQUENTIAL  SHEARS 


A  rotation  of  an  image  through  ancle  ©  implies  that  the  coordinate 
system  is  rotated  through  angle  ©  Then,  the  following  relation¬ 
ship  exists  between  the  nee-  coordinates  (X',  Y')  and  the  old 
coordinates  (X,Y>; 


wh.--e  the  matrix  CR]  indicates  the  resulting  transformation.  If  to  is 
tr a n sf ormat i on  is  divided  into  two  sets  of  operations,  one  along 
scan  lines  and  the  other  along  image  columns,  the  rotation  can  be 
perroraed  in  a  vert  economic.-]  way.  The  "along-line"  operation; 
are  relatively  easy  to  acc omp 1 i sh  The  single  required  input  line 
car,  be  quickly  mapped  to  the  output  image  through  indexing  tables 
computed  from  the  transformation  Efficiencies  in  the  more 
difficult  "a  1  ong -c o  1  umn "  operations  are  achieved  by  moving  groups 
of  columns  as  units,  for  which  multiple  passes  are  made  and  each 
pass  moves  groups  half  as  wide  as  the  last.  This  amounts  to 
performing  the  rotation  by  applying  two  sequential  shearing 
operations,  one  columnwise  and  the  other  rowwise  Sequential 
shearing  directly  achieves  a  satisfactory  rotation  for  small 
angles  However,  for  more  than  a  few  degrees,  the  operation 


introduces  differential  scale  distortions  along  lines  and  columns 
As  with  any  rotation  scheme  applied  to  a  rectangular  pixel  grid, 
resampling  pixels  at  a  compensating  rate  is  required  to  maintain 
scale. 

The  “a  1 ong-c o 1 umn "  operation  is  termed  line  shear  (LS)  and  the 
"along-1 ine"  operation  is  termed  column  shear  <CS)  Then,  the 
transformation  matrix  CR3  is  divided  into  line  shear  and  column 
shear  as  follows: 


nj  Column  she#*  L»n« 


&-cp  a  i 


Of 


The  order  of  line  shear  and  column  shear  can  be  changed  to 
accomplish  the  same  result  as  evident  from  the  following: 
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Ale- or  it  r.ms  io  «- .  c  on  :  1  l  sh  1 1  n  •:  shear  an:  column  shea 
implemented  in  routines  LINil'-n  and  CQLSHR  respectively  A 
interactive  s  o  f  'cmar  e  syste-  r  as  been  implemented  to 
rotation  using  tbe  line  shear  and  column  shear  technique 
software  structure  i:  as  follows: 


r-  ~T ate 


LIhiSHi 


i_SHR 


LINSRR 


3.6.3  LIMITING  ROTATION  ANGLE  USING  SEQUENTIAL  SHEARS 


The  distortion  introduced  by  two  sequential  shears  is  genarally 
negligible  for  small  angles  of  rotation  but  is  quite  significant 
for  angles  which  are  more  than  a  few  degrees,  and  scaling  must  be 
performed  to  restore  image  features  and  their  relationships. 

The  scaling  operation  described  in  the  previous  section  is  very 
nearly  equivalent  to  the  nearest  neighbor  resampling  which  is 
considered  quite  efficient  for  such  operations.  Then#  it  is 
important  to  recognize  that  as  long  as  the  sequential  shearing 
operations  do  not  completely  destroy  the  pixel  neighborhood#  this 
resampling  scheme  has  some  merit. 

During  the  line  shear  and  column  shear  operations#  the  pixels  are 
shifted  around,and  to  limit  their  shifts  to  one  pixel  at  the  most, 
the  rotation  angle  should  be  limited  to  45  degrees  before 
resampling  is  performed.  An  image  of  4096  x  4096  pixels  requires 
12  passes  to  perform  line  shear  of  45  degrees  as  illustated  in 
Figure  3-10. 


3.6.4  ROTATION  ANGLE  OF  90  DECREES 

One  interesting  outcome  of  having  a  general  purpose  line  shear  and 
column  shear  software  capability  is  in  being  able  to  rotate  images 
by  an  angle  of  90  degrees  fairly  easily.  This  is  illustrated  in 
Figure  3-11. 
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Figure  3-10  BINARY  METHOD  OF  PERFORMING  LINE  SHEAR 
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Figure  3-11.  90  Degree  Rotation. 
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Thus,  the  image  is  rotated  by  90  degrees  -rhen  three 
sheering  operations  are  perfc-iiifd  as  folio,  s 
Line  shear  through  4;>  degrees 
Column  shear  through  • ! 5  degrees 
Line  shear  through  degrees 

It  should  be  noticed  that  a  s: cling  operat.-n  is  not 
this  case  and  all  input  dr~a  is  retained  at  the 


sequent : a  1 


needed  in 
end  of  the 


operat i on. 


3.6.5  IMAGE  ROTATION  USING  SEQUENTIAL  SHEARS: 

A  pure  two  dimensional  rotation  about  the  origin  is  performed 
by  the  following  2x2  transformation  matrix: 


1)  when  angle  6  is  positive* 


cos  6 


>m  0 


Sin  6 


Cos.  6 


2)  when  angle  6  is  negative. 


Cos  Q 


&ir\@ 


_Si  r\0 


Cos0 


It  is  evident  that  the  maximum  rotation  required  for  getting 
the  preferred  scan  direction  in  our  mosaicking  approach  is 
190  degrees.  To  accomplish  rotation  using  sequential  shears,  the 
limiting  rotation  angle  requirement  of  section  3.6.3  should  be 
observed  before  resampling  is  performed.  The  following  illust¬ 
rates  the  actual  implementation  of  the  algorithms  to  perform 
image  rotation  to  the  preferred  scan  direction,  image  rectifi¬ 
cation.  and  unrotating  the  corrected  image  to  original  scan 
direction.  The  advantage  of  this  approach  is  that  resampling 
to  be  performed  to  restore  image  features  after  sequential  shears 
is  avoided. 


If  6  is  positive  and  less  than  or  equal  to  45  degrees,  then 
rotation  is  accomplished  by  two  sequential  shears  as  follows: 
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Also<  rotation  in  the  opposite  direction  can  be  implemented  as 
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Nou>  the  image  rectification  process  consists  of 
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where 


£  CS  3  is  column  shear 
C  LS  3  is  line  shear 


The  intent  is  to  have 


C  scaling  ]  [MIC  scaling  ] 
for  R  for  R 

as  one  operation  in  image  correction  software  and  just  perform 
two  sequential  shears  to  rotate  and  un-rotate  image  data. 

This  is  equivalent  to  the  following: 


and  since.  «  0a  »  finally 


This  is  true  when  9  is  negative  and  greater  than  or  equal  to 
-45  degrees.  Thus,  resampling  resulting  from  scaling  is  avoided 

When  0  is  positive  and  greater  that  45  degrees,  then  rotation 
is  accomplished  by  three  sequential  shears  as  follows: 


scaling  column  line  column  rotation 

shear  shear  shear 

step  #4  step  #3  step  #2  step  «1 


Alto,  rotation  in  the  opposite  direction  can  be  implemented  as 


column  line  column  scaling  rotation 

shear  shear  shear 

step  #4  step  #3  step  #2  step  #1 

Then,  the  image  rectification  process  consists  of 

C  CS  3C  LS  DC  CS  DC  scaling  3  C  M  D  C  scaling  DC  CS  DC  LS  DC  CS  D 
with 

C  scaling  ^!  C  M  D  C  scaling  D 
for  R  for  R 

as  one  operation  in  image  correction  software, and  just  perform 
three  sequential  shears  to  rotate  and  un-rotate  image  data. 

This  is  equivalent  to 


and  since.  Q±  =  62  «  finally. 


Thus,  resampling  resulting  from  scaling  is  completely  avoided 


4  THEME  CONTOURING  AND  POLYGON  VECTOR  CONVERSION 


4.1  INTRODUCTION 

Digital  image  file  structures  are  inherently  spatial  in  nature,  and  most  operations 
relating  to  the  extraction  of  spatial  information  are  done  in  a  raster  format  file. 
The  use  of  multiple  data  planes  in  image  files  facilitates  not  only  image  categoriza¬ 
tion,  but  the  interchange  and  logical  manipulation  of  the  resulting  thematic  maps. 
Extracted  information  can  be  rapidly  aggregated  by  various  administrative  and 
physical  units  or  areas  defined  by  map  boundaries.  There  are  other  advantages 
to  raster  format  data  bases,  including  the  ease  of  retrieval  of  spatially 
registered  information,  and  the  ability  to  perform  operations  on  a  pixel  basis 
thus  utilizing  the  finest  resolution  of  the  information  available. 

The  obvious  disadvantage  to  raster  format  data  bases  is  the  great  amount  of 
redundant  storage  used  to  store  smoothly  structured  themes.  This  makes  vector 
format  data  bases,  where  contiguous  areas  are  stored  by  vertex  strings  repre¬ 
senting  outlines,  very  attractive.  Relative  storage  economy  for  each  format  is 
determined  by  the  ratio  between  the  effective  area  subtended  by  a  pixel  and  the 
average  size  of  uniformly  mapped  areas.  The  advantage  of  vector  storage  grows 
exponentially  with  area  size,  showing  an  advantage  over  raster  storage  when  the 
average  contiguous  area  is  about  twenty  five  pixels.  Vector  format  data  bases 
are  generally  more  economical,  and  are  particularly  easy  to  use  when  generating 
chloroplethic  maps  from  stored  information.  The  problem  is  that  the  original 
extraction  of  information  to  be  mapped  is  very  difficult. 

The  Standard  Interface  Format  (SIF)  provides  a  common  mechanism  for  the  transmittal 
of  graphics  and  spatial  information  between  dissimilar  data  bases.  Information  is 
passed  in  vector  format  on  a  magnetic  tape  containing  one  or  more  files.  Each  file 
is  used  to  describe  one  drawing  or  theme. 


SIF  FILt  FORMAT 


An  SIF  tape  is  made  up  of  one  or  more  files.  Each  file  represents  one  drawing 
(or,  for  our  purposes,  one  theme). 

Each  file  contains  one  or  more  data  blocks,  and  each  data  block  is  256  32-bit 
words  in  length.  The  first  word  of  each  block  is  always  a  command  header. 

Data  follow  the  header. 

COMMAND  AND  DATA  FORMAT  FOR  LINE  STRINGS 

The  first  word  of  each  command  is  the  command  header.  It  contains  four  bytes 
reflecting  the  type  of  command,  its  size,  a  command  modifier,  and  a  one  byte 
value  in  some  cases.  The  following  is  the  SIF  format  for  line  strings: 


COM1AND  -  40  for  a  line  string 

-  80  for  continuation  of  a  line  string 

#OF  WORDS  -  Number  of  32-bit  words  in  the  command.  (255  maximum.) 

SUBTYPE  -  1  for  open  line  string 

-  2  for  closed  line  string  representing  an  area 

-  3  for  closed  line  stirng  representing  a  hole 

VALUE  -  0  for  line  strings 
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The  mechanics  of  converting  themes  to  polygons  involves  an  unknown  and  potentially 
large  number  of  vertices  and  presents  serious  problems  concerning  computer  storage 
requirements.  The  number  of  vertices  needed  to  outline  classes  in  a  typical  512  x 
512  pixel  scene  can  be  grossly  estimated  at  between  5000  and  15000,  with  the  number 
being  highly  sensitive  to  overall  scene  complexity.  In  an  operational  situation 
where  an  average  map  sized  area  is  being  processed  at  high  resolution,  larger 
images  or  image  mosaics  of  4000x  4000  pixels  may  be  a  more  representative  processing 
unit  size.  We  clearly  need  efficient  vertex  storage,  methods  of  recycling  storage 
during  conversion,  and  a  means  of  dissecting  large  files  into  smaller  independent 
segments  for  processing.  Directly  storing  vertex  coordinates  in  SIF  values  would 
use  sixty-four  bits  per  vertex,  and,  if  workable,  would  be  very  extravagant. 

Chain  codes  are  an  exceptionally  efficient  method  for  storing  vertex  lists,  particul 
ly  if  the  polygons  are  of  complex  shapes  containing  only  short  straight  line  seg¬ 
ments.  A  chain  code  is  comprised  of  a  sequence  of  single  digit  numbers  0-7,  one 
number  per  vertex.  Based  on  a  rectangular  pixel  grid,  each  of  the  numbers  defines 
one  of  the  eight  possible  directions  pointing  to  the  next  vertex.  The  complete 
geometric  description  of  any  polygon  is  provided  by  the  coordinates  of  some  start¬ 
ing  point  and  the  polygon's  chain  code.  There  is  some  waste  associated  with  storing 
redundant  vertices  along  extended  straight  line  segments,  but  this  is  unimportant 
in  view  of  the  economy  gained  in  storing  an  average  polygon.  With  the  restricted 
set  of  direction  code  possibilities  only  four  bits  are  required  to  store  each 
element  in  the  chain.  In  this  respect  we  might  compromise,  wasting  four  bits  per 
vertex  by  using  bytes  to  decrease  execution  time  and  program  complexity. 

Careful  consideration  must  also  be  given  to  the  definition  of  "outline."  When 
dealing  with  a  continuous  surface  in  three  space,  the  concept  of  a  contour  is 
straightforward  and  can  be  defined  as  the  intersection  of  that  surface  with  some 


plane  within  the  space.  In  the  discrete  world  of  the  pixel  grid  used  with  raster 
format  data  bases,  the  concept  of  a  contour  is  a  little  more  ambiguous.  The  idea 
of  an  intersecting  plane  is  insufficient  as  the  plane  no  longer  intersects  a 
continuum.  Actually ,  when  dealing  with  pixels,  a  plane  is  no  longer  considered 

t 

continuous,  but  a  series  of  discrete  locations  spaced  at  regular  intervals.  Any 
surface  intersection  in  a  discrete  space  must  also  be  discontinuous,  as  pixel 
positions  represent  the  limit  of  available  resolution. 

Several  definitions  pertaining  to  discrete  geometry  are  in  order: 

DIRECT  NEIGHBORS  -  Pixels  pairs  that  are  directly  adjacent,  sharing  a  common 
side  in  the  pixel  grid. 

INDIRECT  NEIGHBORS  -  Pixel  pairs  that  share  a  common  corner  point  in  the  pixel  grid. 

(N)-NEIGHBOR  -  The  pixel  neighbor  in  direction  N  from  any  given  position.  N 
ranges  from  0  to  7  and  increases  with  counter  clockwise  angle  from  east. 

DIRECT  PATH  and  INDIRECT  PATH  -  Sequences  of  pixels  containing  direct  neighbors, 
and  any  neighbors,  respectively. 

INDIRECTLY  CONNECTED  REGION  -  A  set  of  pixels,  each  pair  of  which  can  be  connected 
by  a  path  of  direct  or  indirect  neighbors.  Similarly  a  DIRECTLY  CONNECTED  REGION 
is  a  set,  each  pair  of  which  is  connected  by  a  path  of  direct  neighbors. 

CONTOUR  or  INDIRECT  CONTOUR  -  The  set  of  all  pixels  in  a  region  which  have  at 
least  one  direct  neighbor  not  in  the  region, 

DIRECT  CONTOUR  -  The  set  of  all  pixels  in  a  region  which  have  at  least  one  neighbor 
not  in  the  region. 

As  the  list  indicates,  there  are  several  possible  contour  definitions  usable  for 
theme  outlining,  including  direct  and  indirect  paths  outlining  direct  and  indirectly 
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connected  regions.  Our  choice  is  dictated  by  data  storage  economy.  We  need  to 
describe  themes  with  the  fewest  number  of  vertex  strings,  and  to  use  the  shortest 
strings  possible  to  provide  a  complete  description.  Of  the  possible  methods, 
the  one  best  suited  for  outlining  noisy  or  irregular  themes  is  to  describe 
indirectly  connected  theme  areas  with  indirect  contours. 

The  method  for  outlining  theme  holes  must  also  be  established.  If  we  are  willing 
to  have  hole  outlines  occur  within  the  set  defined  by  a  theme  region,  then  the 
same  contour  definition  can  be  used  for  holes.  This  keeps  the  software  simpler. 

A  single  contouring  algorithm  utilizing  an  operator  which  marches  along  the  theme 
boundary  and  selects  the  rightmost  available  next  pixel,  can  be  used  for  both 
region  types.  The  types  are  even  classified  by  the  direction  of  outline  traverse. 
The  single  algorithm  scheme  does  create  some  special  and  potentially  troublesome 
conditions  which  must  be  remembered.  First,  since  the  hole  contours  will  occur 
on  theme  pixels,  calculation  of  hole  areas  must  be  different  from  theme  areas  to 
compensate  for  peripheral  difference  of  one  pixel.  Also,  contours  of  the  two 
region  types  have  different  spatial  characteristics.  The  theme  outlines  can 
contain  degenerate  polygon  segments  which  have  no  width,  and  hence  no  interior 
area.  Examples  are  single  pixel  width  protrusions,  single  pixel  width  bridges 
between  larger  areas,  and  also  isolated  lines  and  single  pixels.  All  of  these 
features  can  be  considered  degenerate,  spanning  no  interior  area  in  the  discrete 
plane.  Hole  contours,  if  they  are  defined  to  be  in  the  theme  set,  are  a  little 
easier  to  work  with  as  they  cannot  contain  degenerate  segments.  Another  commonly 
occurring  special  condition  arises  when  a  theme  area  includes  a  hole  that  is 
positioned  within  one  pixel  of  the  region  boundary.  In  these  cases  the  associated 
contours  will  be  coincident  for  some  fraction  of  their  length,  and  can  complicate 
processing. 


4.2  THEME  OUTLINING  ALGORITHM 


The  outlining  algorithm  is  a  type  of  marching  operator,  and  as  such  must  be 
capable  of  handling  any  outline  or  hole,  providing  a  predictable  and  repeatable 
outline  for  any  input  theme  of  any  complexity.  With  the  algorithm  being  used 
hundreds  or  even  thousands  of  times  on  a  single  complex  image,  an  operator  which 
occasionally  gets  off  track  or  fails  to  close  outlines  in  complex  situations  is 
unacceptable.  Fortunately,  such  algorithms  have  been  developed  and  have  evolved 
to  a  mature  state.  The  one  chosen  is  a  contouring  algorithm  which  basically 
starts  at  some  point  on  a  region  contour,  determines  a  starting  direction  and 
marches  along,  finding  the  rightmost  next  pixel  satisfying  the  contour  definition. 

The  actual  pixel  testing  sequence  is  interesting  and  intricate.  Directions  close 
to  the  previous  one  are  the  most  probable  for  contour  continuation  and  are  tested 
first.  Leftward,  but  not  rightward , acute  angle  turns  are  allowed.  This  results 
in  the  generation  of  an  indirect  contour  around  theme  areas.  Figure  4-1  is  a 
pseudo  code  algorithm  description.  Some  complicating  factors  are  included  to 
handle  discrete  plane  contours.  A  loop  counter  is  needed  to  count  the  right 
turns  made  at  any  particular  operator  position  in  order  to  identify  single,  isolated 
pixels  and  to  avoid  continually  circling  on  them.  If  no  new  contour  pixel  is 
found,  the  tests  used  are  such  that  the  marching  operator  will  retrace  its  steps 
until  a  new  rightward  pixel  is  found.  This  feature  makes  it  possible  to  handle 
single-pixel-width  lines  or  outline  segments.  Finally,  a  dual  test  is  used  for 
contour  termination.  The  operator  must  return  to  the  start  point  and  must  be 
heading  in  the  original  starting  direction.  The  reason  behind  the  first  test  is 
obvious,  and  the  second  is  to  avoid  premature  termination  when  the  starting  point 
happens  to  lie  on  a  section  of  the  polygon  where  the  operator  is  retracing  its  steps. 
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The  outline  of  *  connected  region 
1*  completely  specified  by: 

Starting  Pixel  -  <x,y)  ,and 
Chain  Code  -  5757802131 15543 


The  Outline  of  a  connected  region  R  is  defined  as  the  set  of  pixels  in 
R  shieh  hsve  at  least  one  direct  neighbor  (one  sharing  a  coascn  side) 
not  in  R. 


to  outlining  operator  can  be  described  as  sose thing  that  steps  along  pixel 
to  pixel  on  the  edge  of  a  region  by  selecting  the  rightmost  available 
satisfying  the  outline  definition. 


Theses  are  considered  to  have  a  positive  area#  and  are  outlined  in  a  counter 
clockwise  direction.  Hales  are  outlined  eloekoioe. 


Defining  an  outline  oith  a  chain 
Figure  4-1. 


The  contour-tracing  algorithm  is  presented  in  Figure  4-2.  However,  the  completely 
detailed  algorithm  is  given  in  the  Fortran  program  listing  for  the  "TRACKR" 
module  which  is  largely  self-documenting.  For  input,  it  requires  only  that  a 
start  point  be  selected  which  is  within  the  region  being  contoured  while  its 
directly  adjacent  neighbor  to  the  west  is  not.  Output  of  the  module  is  a  closed 
contour  chain  code  loaded  into  a  queue.  An  image  I/O  subroutine,  "UPDATE"  is 
called  by  TRACKR  and  does  several  things.  It  accepts  the  just  located  direction 
to  the  next  contour  pixel  and  reads  the  adjacent  image  line  in  that  direction, 
using  a  rotating  buffer  scheme  for  efficiency.  It  returns  values  for  the  eight 
pixels  surrounding  the  next  operator  position.  It  also  marks  the  pixel  just 
departed  by  incrementing  its  value  by  one,  thus  permanently  identifying  the 
contour  as  having  been  traced.  Troublesome  single-pixel -width  segments  are 
uniquely  marked,  as  their  values  will  be  incremented  twice.  In  an  effort  to 
minimize  the  large  number  of  image  line  reads  and  writes,  a  buffering  scheme  is 
used  which  holds  the  three  lines  surrounding  the  operator,  reading  or  writing 
single  lines  when  the  operator  moves  up  or  down  on  the  image. 

At  this  point  it  is  appropriate  to  refer  to  Figure  4-3  which  presents  all  the  con¬ 
touring  and  polygon  conversion  modules  and  their  inter-relationships. 

SIFDMP  extracts  numeric  contour  coordinates.  It  is  an  optional  driver  for  the  sub¬ 
routine  TRACKR,  with  expanded  output  capabilities,  and  is  used  to  generate  a 
single  theme  SIF  magnetic  tape.  SIFDMP  is  called  once  for  each  region  in  the  theme. 
Given  a  contour  start  point  in  the  same  manner  as  TRACKR,  this  subroutine  similary 
locates  the  vertices,  counts  them  to  insure  the  SIF  line  string  minimum  of  three 
is  surpassed,  then  converts  them  into  strings  of  thirty-two  bit  cartesian  coordinates 
in  accordance  with  SIF  record  structure.  Record  header  information  is  assigned 
to  prefix  the  vertex  strings  and  each  record  is  written  to  tape  after  it  becomes 
filled.  An  added  feature  that  was  included  but  not  currently  used  is  the  capab¬ 
ility  to  linearly  transform  the  polygon  vertices  into  any  arbitrary  coordinate 
system. 

4-8 


STARTING  AT  SOM  POINT  A  IN  REGION  R  8UX  THAT  ITS  4-fCIMOR  IS  NOT  IN  R 

SET  TX  OARENT  PI*L  TO  T>C  INITIAL  PDG. 

C-A 

FIRST  ■  .TALE. 

SET  INITIAL  SEARCH  DIRECTION  S-6 

WHILE  C(C  .ft.  A)  .AM).  (SLOPE  .ft.  D)>  .OR.  (FIRST  •  .TRUE.)  DO 
FXUC)  -  .FALSE. 

OOINTER  -  1 

WILE  (.NOT.  FDUM»  .AM).  (COINTER  .LE.  3)  DO 

IF  (B,  TVE  (S-l)-tEIGVDOR  OF  C  IS  IN  R) 

TVEN: 

nut  •  .me. 

C-B 

S-S-2 

ELSE  IF  (B,  TM  (S)-fCIMOR  OF  C  IS  IN  R) 

F(Ut  •  .TRUE. 

C-B 

OX  IF  (B,  TVE  (SfX)-fCI©*0R  OF  C  IS  IN  R) 

FCLND  -  .TRUE. 

C-B 

ELX:  INCRDCNT  SEARCH  DIRECTION  BY  SB  KGREES 
S-S+2 

COINTER  •  COINTER  41 
DC  IF 

EM)  WILE 

ADD  DIRECTION  CODE  OF  LOCATED  PDCL  TO  CHAIN 

IF  (FIRST  -  .TRUE.) 

TMN:  D  -  DIRECTION  CODE  FOR  B 
OX:  SLOPE  -  DIRECTION  CODE  FOR  B 
EM)  IF 

END  WILE 


A  -  Starting  point  on  the  contour  of  region  R 

C  -  Curront  point  Woo*  neighborhood  i*  being  ex  Mined 

S  -  Pixel  neighboring  C 

D  -  Slope  of  contour  leaving  point  A 

ftlthortlc  it  oodulo  8 


Contour  Treeing  Algor itho 

Figure  4-2. 


TFWCXR 


CNTOUR  -  M»ln  program  xhlch  provides  an  operator  Interface  and  controls  the 
selected  processing. 

REGION  -  Defines  the  scene  segment  to  he  processed  by  Inserting  a  rectangular 
window  Into  the  scene.  The  window  will  be  searched  Internally  for 
the  location  of  theme  areas. 

TRACW?  -  The  contour  tracing  subroutine.  Given  some  point  on  the  outline 
of  a  theme  area  or  hole,  this  subroutine  traces  an  outline  and 
stores  the  associated  chain  code. 

SOAMER  .  Region  scanning  routine  a#ilch  systematically  searches  the  Interior 
of  a  designated  these  area  or  hole,  and  returns  the  location  of 
any  internal  holes  or  areas.  The  Interiors  of  any  features  located 
are  not  scanned. 

SIFDTP  -  A  subroutine  ah ich  generates  Sir  tape  records  from  a  region  chain 
code  generated  by  TRAO0?. 

UPDATE  *  Handles  1  stage  file  I/O  required  by  the  contour  tracing  subroutines. 
After  each  new  pixel  position  Is  found  by  90ATCR,  this  routine 
marks  Its  position  and  returns  the  surrounding  pixels. 

IRK  -  G.E.  DIAL  cursor  position  reading  subroutine. 

IRV  -  G.E.  DIAL  I  stage  line  reading  subroutine. 

IWV  -  G.E.  DIAL  image  line  writing  subroutine. 


Theme  Contouring  and  Polygon  Conversion  Program  Nodules 

Figure  4-3. 


4.3  THEME  REGION  LOCATING  ALGORITHM 

The  vector  conversion  software  package  has  a  second  major  module  used  to  locate 
and  tabulate  theme  areas.  It  is  a  region  scanner  which  takes  any  region  in  a 
thematic  map  and  systematically  searches  the  interior  for  the  presence  of  sub- 
regions.  The  subregions  may  be  either  theme  areas  or  holes,  depending  on  the 
character  of  the  region  itself.  But  in  order  to  keep  track  of  region  nesting, 
it  is  important  that  the  search  area  is  internal  to  the  region  and  external  to 
any  located  subregions.  Generations  of  subregions  then  can  be  tracked  by  the 
number  of  times  the  scanning  algorithm  must  be  called  to  reach  a  given  level  of 
nesting. 

Coupled  with  the  previously  described  outlining  algorithm,  the  region  scanner 
provides  a  means  of  locating  and  extracting  all  contours  associated  with  a  single 
theme.  If  desired,  a  tree  structure  could  be  used  to  keep  track  of  which  holes 
are  contained  in  which  areas,  etc.  Since  the  outlining  algorithm  easily  re¬ 
generates  vertex  strings  completely  describing  a  region  from  its  starting  point, 
the  data  structure  referencing  all  theme  regions  need  only  store  a  single  point 
for  each  region. 

Figure  4-4  shows  the  scanning  algorithm  and  sequence  used  in  searching  for  holes 
in  a  theme  area.  A  slightly  simplified  hole-finding  algorithm  is  presented  for 
clarity.  The  generalized  algorithm  is  found  in  the  Fortran  listing  for  SCANER. 
The  algorithm  requires  some  point  on  the  outline  of  some  region,  draws  the 
complete  outline,  and  then  finds  all  areas  of  the  opposite  parity  (areas/holes) 
enclosed  in  the  region.  The  module  returns  start  points  for  all  subregions  in 
direct  contact  with  the  region  being  scanned. 
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GIVEN  A  POINT  ON  TVE  BOUNDARY  OF  7>C  REGION  TO  BE  INTERNALLY  9CANCD 
GftLL  TRACXR  TO  NARK  TVE  REGION  OUTLIN  AN)  STORE  IT  IN  ft  QUEUE 

WILE  (QUEUE  NOT  EMPTY)  DO 

SET  LAST  ODIN  CODE  TO  CURRENT  C0-C 

GET  DIRECTION  OPPOSITE  TO  C0  ftCB-NOO( (C0*4),B) 

GET  NXT  CHAIN  CODE  FROM  QUEUE  C -QUEUE (FIRST  ELEMENT) 

TEST  FOR  DESCENDING  PftTH  BEFORE  SCATHING  LEFT  TO  RIGHT 
IF  ((C.N.0).AN).(C8.N.4).AN>.(AC8.GE.C)>  TVEN: 

TERMINATE  •  .FALSE. 

WILE  (SCAN  NOT  TERMINATED)  DO 

EXTRACT  NEXT  SEQUENTIAL  SET  OF  TVREE  PI«LS  (ABC)  FROM  SCAN  LINE 

TEST  FOR  INTERSECTION  WITH  ft  PREVIOUSLY  IN.  OCA  TED  SUBREGION 
IF((A-0).AND. (B-l) )  .OR.  ((A-0).AN>.(B-2).AN>.(C-0>) 

TVEN: 

SAVE  NOTOLK)  SUBREGION  START  POINT  FOR  FUTIRE  REFERENCE 
CALL  TRAOR  TO  MARK  SUBREGION  AN)  APPEND  IT  TO  QUEUE 
TEW1INATE  •.TRUE. 

ELSE: 

TEST  FOR  ft  PREVIOUSLY  MAR>0  OUTLIfC 
IF  (B.GE.2)  TEJWINATE  -.TRUE. 

EN>  IF 
EN)  WILE 
EN)  IF 
EN)  WILE 


—Region  Scanning  Algorithm  To  Locate  Hole* — 

Figure  4-4. 
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The  procedure  starts  with  a  binary  image  in  which  theme  pixels  have  been  assigned 
a  value  of  one,  and  "not  theme"  pixels  are  zero.  The  outline  start  point  is  used 
to  call  the  contour  tracing  module,  which  marks  each  pixel  on  the  outer  boundary 
by  adding  one  to  the  pixel  values.  Since  the  outline  traced  is  wholly  included 
in  the  theme,  outline  pixels  will  now  have  a  value  of  two  and  thus  be  distinguished 
from  untraced  outline  pixels.  As  the  outline  is  traced,  the  resulting  chain  code 
is  temporarily  placed  in  a  large  queue  in  preparation  for  scanning  the  interior 
of  the  region. 

Points  for  starting  interior  scans  are  found  in  the  queue.  Chain  code  elements 
are  extracted  one  by  one,  and  sequential  pairs  are  tested  to  determine  if  they 
lie  on  ascending  or  decending  segments  of  the  region  outline.  Ascending /descending 
path  pixels  are  used  as  start  points  for  rightward  scans  searching  for  areas/holes. 
Along  the  scan,  location  of  either  type  subregion  is  found  by  the  occurrence  of 
two  pixels  having  a  value  sequence  of  zero,  one.  This  sequence  occurs  where  the 
scan  line  intersects  the  leftmost  edge  of  a  theme  (climbing  on),  or  the  rightmost 
edge  of  a  hole  (climbing  out).  An  additional  sequence  of  zero,  two,  zero  locates 
holes  whose  right  side  is  coincident  with  the  right  side  of  the  region  being 
scanned.  Immediately  on  location  of  a  subregion,  the  current  scan  coordinates 
are  stored  for  future  reference. 

The  contour  tracer  is  called  to  mark  the  subregion  outline,  and  to  append  the 
outline  to  the  queue.  A  particular  scan  is  stopped  when  a  new  subregion  is 
located,  when  the  outline  of  previously  traced  subregion  is  found,  or  when  the 
traced  right  side  of  the  region  is  reached.  Both  the  latter  cases  are  identified 
by  pixel  values  of  two  or  greater,  as  they  will  have  been  traced  and  marked  at 
least  once. 


The  reason  for  using  a  queue  to  store  region  chain  codes  now  becomes  apparent. 

The  region,  if  it  is  a  theme,  will  be  outlined  in  a  counter-clockwise  direction. 
Scanning  from  descending  segments  will  locate  holes  to  the  right  (provided  they 
are  not  eclipsed  by  other  holes)  and  immediately  enter  them  into  the  queue. 
Eventually,  when  the  region  start  point  is  returned  to,  the  region  outline  will 
have  been  completely  removed  from  the  queue.  The  scanning  algorithm  continues, 
with  the  upcoming  queue  elements  belonging  to  the  outline  of  the  first  hole. 

Since  holes  are  outlined  in  a  clockwise  manner,  descending  segments  occur  on 
their  right  sides,  conveniently  serving  as  start  points  for  scans  into  sections 
of  the  region  shaded  from  the  left  by  the  holes.  By  the  time  the  queue  is  finally 
empty,  the  region  boundary  along  with  all  hole  outlines  will  have  been  checked 
for  descending  segments.  Consequently  all  sections  of  the  theme  region  external  to 
the  hole  type  subregions  will  have  been  scanned. 

The  scanning  algorithm  presented  in  Figure  4-4  must  be  expanded  somewhat  to  find 
areas  within  holes  in  order  to  become  a  generalized  region  scanning  tool.  The 
decision  as  to  where  on  the  chain  to  start  a  line  scan  becomes  more  complex.  To 
preserve  the  left-to-right  scanning  convention  when  searching  for  theme  areas, 
we  modify  the  logic  to  scan  from  ascending  chain  segments  (as  opposed  to  the 
descending  segments  used  in  hole  location).  Figure  4-5  shows  the  derivation  of 
a  comprehensive  scan  start  test,  useful  for  identifying  either  ascending  or  decend- 
ing  segments.  The  test  also  locates  chain  code  turns  at  local  tops  and  bottoms 
in  outlines  where  a  simple  ascent/descent  test  is  insufficient  to  identify  start 
points.  In  all  cases  the  location  of  a  previously  traced  outline  is  used  to 
terminate  the  scan. 


Starting  scans  from  ascending 
segments  of  these  outlines 


CB  cases 


Corresponding  C's 
for  scan  start 


Notation:  C  -  Current  chain  code  pointing  to  next  position 

C0  -  Last  chain  code  pointing  to  current  position 
ROB  -  Direction  opposite  to  CB 


Rules: 


1.  Never  scan  «4>en  CB*B  as  the  line  has  already  been  scanned 

2.  Never  scan  «ten  Ca4  as  a  better  start  position  is  coming  up. 

3.  For  rightward  scan  frost  ascending  chain  segments,  C  oust  be 
included  in  a  counter  clockwise  arc  between  fC  (1)  and  RC8. 

4.  For  rightward  scan  frost  descending  chain  segswnts,  C  oust  be 
included  in  a  clockwise  arc  between  SE  (7)  and  ROB. 


Ascending  segment  scan  start  test  —  (C.NE.B).ftQ. (CB.fC.4).flfC. (AC0.GE.C) 
Descending  segment  scan  start  test  —  (C.fC.B).flN). (CB.fC.4).PTO. (RCB.LT.C) 


Conditions  for  Starting  External  Region  Scans 

Figure  4-5. 


A  final  complication  concerns  the  handling  of  degenerate  outline  segments.  When 
searching  for  theme  are  ,  scanning  is  internal  to  a  hole  type  region,  hence 
external  to  any  located  heme-type  subregions.  If  the  region  is  to  be  completely 
scanned,  all  ascending  egments  including  degenerate  ones  must  be  used  as  scan 
start  points.  For  noisy  themes  it  is  particularly  important  to  scan  from  singeltons 
that  can  otherwise  have  a  combining  effect  that  eclipses  large  areas.  Fortunately  the 
exterior  region,  being  a  hole,  will  not  contain  degenerate  sections  and  thus 
allows  the  entire  queue  contents  to  be  processed  consistently.  For  the  alternate 
case,  the  search  is  internal  to  a  theme  area  and  scanning  from  degenerate  segments 
must  be  suppressed  or  the  scan  will  be  outside  of  the  region.  Figure  4-6  shows 
the  possibilities. 

The  scanning  algorithm  has  one  known  limitation.  If  SCANER  is  being  used  to  locate 
holes  (a  second  level  operation  for  our  purposes),  and  if  a  hole  contains  a  nested 
theme  subregion,  and  if  all  portions  of  that  hole  on  image  lines  above  the  nested 
theme  are  scan  shaded  from  the  left  by  other  holes,  an  irregularity  occurs.  As 
the  lines  being  scanned  descend  past  the  bottom  of  the  shaded  zone,  the  next  scan 
will  enter  the  previously  untraced  hole,  and  instead  of  finding  the  far  side,  will 
prematurely  trace  around  the  nested  theme.  Figure  4-6b  illustrates  the  error  pro¬ 
ducing  theme  configuration.  Tracing  direction  will  still  be  correct,  but  without 
more  intricate  logic  the  nested  theme  will  be  associated  with  the  wrong  external 
region.  For  our  immediate  thematic  mapping  application  this  irregularity  will 
probably  go  unexercised  as  it  can  not  occur  without  a  rare  three-level  nesting. 

The  remaining  module  "REGION"  is  used  to  insert  a  hole  into  the  stored  image  to 
define  the  processing  window.  A  start  point  on  the  window  is  returned  and  sub¬ 
sequently  passed  to  SCANER  to  begin  processing  of  the  enclosed  theme  areas. 


Level  1  region  scan  from  ascending 
queue  segments  to  locate  theme 
areas. 


Scans  are  internal  to  hole  region 
and  external  to  located  areas. 


•  Start  scans  from  single  pixel 
width  segments.* 


•  Start  scans  from  singletons. 


Level  2  scan  of  decending  queue 
segments  to  locate  holes  within 
theme  areas. 


Scans  are  internal  to  theme  area 
and  external  to  located  holes. 


•  Don't  scan  from  single  pixel 
width  segments 


•  Don't  scan  from  singletons. 


Modifications  of  left-to-right  scan  start  logic  for  handling  degenerate  polygons. 


Figure  4-6. 


(X,  -  ax)  X,  +  (V,  -  ay)  y, 


(X2  -  ax)  x2  +  (Y2  -  ay)  y2 


Cross-multiplying,  we  get: 

{(V,  -  ay)  X,  -  (X,  -  ax)  y,}  fa  -  ax)  x2  *  (V2  -  ay)  y2j  = 

{(X,  -  ax)  x,  +  (Y,  -  ay)  y,j  ^  -  ay)  x2  -  (X2  -  ax)  y2^ 

Expanding, 

(Y1  -  ay)  (X2  -  ax)  x1  x2  +  (Y1  -  ay)  (Y2  -  ay)  x1  y2  -  (X1  -  ax)  (X2  -  ax)  x. 
-  (X1  -  ax)  (Y2  -  ay)  y1  y2  =  (Y2  -  ay)  (X1  -  ax)  x1  x2  + 

(Y2  -  ay)  (Y,  -  ay)  X2  y,  -  (X2  -  ax)  (X,  -  ax)  y2  x,  -  (X2  -  ax)  (Y,  -  ay) 

Collecting  terms  in  x1  x2,  x1  y2,  x2  y^,  and  y1  y2  together 

{(X,  -  ay)  (X2  -  ax)  -  (Y2  -  ay)  (X,  -  a^x,  x2  + 

{(Y,  -  ay)  (Y2  -  ay)  +  (X2  -  ax)  (X,  -  a^x,  y2  - 

{(X1  '  ax*  <X2  ’  ax*  *  (Y1  '  V  (V2  "  ay>}  x2  y1  ‘ 

/(X,  -  a,)  (Y,  -  a  )  -  (X.  -  a,)  (Y,  -  av)}y,  y,  =  0 


{(Y1  -  ay)  (X2  -  ax)  -  (Y2  -  ay)  (X1  -  ax)}  (x1  x2  +  y,  y2)  + 

^(Y1  -  ay)  (Y2  -  ay)  +  (Xj  -  ax)  (X2  -  ax)J  (x1  y2  -  x2  y^  =  0 

Dividing  by  (x^  y2  -  x2  y^)  throughout,  we  get 

\  X1  x2  "  y1  y2 

VY1  "  ay)  ~  ax}  ”  ^Y2  ‘  ay^  ^X1  ‘  ax$  ^  X1  y2  -x2  y1  ^ 

{(V,  -  ay>  <y2  -  ay)  +  (X,  -  ax)  (X2  -  ax)]  =  0 

provided  x1  y2  t  x2  yr  the  contro1  points  CPI  and  CP2  and  the  origin 

y1  y2 

of  the  photo-coordinate  system  are  non-col  inear,  then  t  ~  and 

X1  x2 

X,  y2  -  x2  y,  4  o. 

Let  us  call 

(x1  x2  +  y1  y2)/(x1  y2  -  x2  y^  =« 

Then, 

(Y1  X2  -  ax  V1  -  ay  h  +  a*  ay  -  V2  X1  +  ax  V2  +  ay  *1  ‘  ax  ay}~ 

+  (x,  X2  -  ax  (X,  +  X2)  +  ax'  +  Y,  Y2  -  ay  (Y,  t  Y2)  *  ay’}  =  0. 

Collecting  terms  in  a  and  a  , 

X  Jr 

ax*  +  ay'  -  ax  {(Y,  -  Y2)~  +  (X,  +  X2^  -  ay  ((X2  -  X, )«  +  (Y,  +  y2)| 

+  {x,  x2  +  y1y2+(y,x2-y2x1)^  =  0. 


A-2 


or 


V  ♦  V  4  ax  (<v2  '  V*  -  (X1  *  X2>}  4  ay 

4  fx1X2  +  Y1Y2t(Y1X2-Y2X1)~}=  0. 

or 

ax*  +  ay2  +  JTax  +  *?ay  +  £  =  0. 

where 

*  =  (Y2  -  -  (X1  +  x2) 

^  =  (X1  -  X2)o<  -  (Y1  +  Y2) 

*■  Xl  X2  4  V1  V2  +  <XI  X2  '  h  Xl>« 


A-3 


